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We study, using numerical simulations, the dynamical evolution of self-gravitating point particles 
in static euclidean space, starting from a simple class of infinite "shuffied lattice" initial conditions. 
These are obtained by applying independently to each particle on an infinite perfect lattice a small 
random displacement, and are characterized by a power spectrum (structure factor) of density 
fluctuations which is quadratic in the wave number fc, at small k. For a specified form of the 
probability distribution function of the "shuffling" applied to each particle, and zero initial velocities, 
these initial configurations are characterized by a single relevant parameter: the variance 5^ of the 
"shuffling" normalized in units of the lattice spacing I. The clustering, which develops in time 
starting from scales around I, is qualitatively very similar to that seen in cosmological simulations, 
which begin from lattices with applied correlated displacements and incorporate an expanding spatial 
background. From very soon after the formation of the first non-linear structures, a spatio-temporal 
scaling relation describes well the evolution of the two-point correlations. At larger times the 
dynamics of these correlations converges to what is termed "self-similar" evolution in cosmology, 
in which the time dependence in the scaling relation is specified entirely by that of the linearized 
fluid theory. Comparing simulations with different 5, different resolution, but identical large scale 
ffuctuations, we are able to identify and study features of the dynamics of the system in the transient 
phase leading to this behavior. In this phase, the discrete nature of the system explicitly plays an 
essential role. 
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I. INTRODUCTION 



The problem of the evolution of self-gravitating clas- 
sical particles, initially distributed very uniformly in in- 
finite space, is as old as Newton. Modern cosmology 
poses essentially the same problem as the matter in the 
universe is now believed to consist predominantly of al- 
most purely self-gravitating particles — so called dark 
matter — which is, at early times, indeed very close to 
uniformly distributed in the universe, and at densities 
at which quantum effects are completely negligible. De- 
spite the age of the problem and the impressive advances 
of modern cosmology in recent years, our understanding 
of it remains, however, very incomplete. In its essen- 
tials, i.e., stripped of the full detail of current cosmolog- 
ical models, it is a simple well posed problem of out of 



equilibrium statistical mechanics^. In this context, how- 
ever, it has been relatively neglected, primarily because 
of the intrinsic difficulties associated with the attractive 
long-range nature of gravity and its singular behavior at 
vanishing separation. In recent years there has, how- 
ever, been renewed interest (see e.g. [^) in the physics 
of systems with long-range interactions, in which context 
self-gravitating systems are one of the paradigmatic ex- 



Strictly speaking it is not actually known whether the problem 
is well-controlled without a regularizatiorLof the singularity in 
the gravitational force at r = (see e.g. for a recent discus- 
sion and list of references). In practice, in numerical simulation, 
there is no intrinsic problem in implementing the N-body gravi- 
tational dynamics without such a regularisation for typical initial 
conditions (i.e. in which particles are not placed initially at the 
same point). In the numerical simulations reported here, as in 
cosmological simulations, we do, however, use such a regulariza- 
tion. This is done solely for numerical efficiency, and the results 
analysed are tested numerically for their independence of the 
associated cut-oflf (see below). 
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amples (for a review, see e.g. |^). A considerable amount 
of work on these systems in this context has focussed on 
finite systems (see e.g. [§, |], ^, 0|) — in which, in certain 
cases, some of the instruments of equihbrium statistical 
mechanics may be applied ^ — and on more tractable 
onc-dimcnsional models (see e.g. ^ ^). 

In cosmology perturbative approaches to the problem, 
which treat the very limited range of low to modest am- 
plitude deviations from uniformity, have been developed 
(see e.g. ^^), but numerical simulations are essen- 
tially the only instrument beyond this regime. While 
such simulations constitute a very powerful and essential 
tool, they lack the valuable guidance which a fuller ana- 
lytic understanding of the problem would provide. The 
dynamics of infinite self-gravitating systems is thus both 
a fascinating theoretical problem of out of equilibrium 
statistical mechanics, directly relevant both in the con- 
text of cosmology and, more generally, in the physics of 
systems with long-range interactions. 

Approaching the problem in the context of statistical 
mechanics, as we do here, it is natural to start by reduc- 
ing as much as possible the complexity of the analagous 
cosmological problem. We wish to focus on the essen- 
tial aspects of the problem. Thus we consider clustering 
without the expansion of the universe, and starting from 
particularly simple initial conditions. With respect to the 
motivation from cosmology, there is of course a risk : in 
simplifying we may loose some essential elements which 
change the nature of gravitational clustering. Our results 
suggest that this is not the case. Even it were, it seems 
unlikely that we will not learn something about the more 
complex cosmological problem in addressing this slightly 
different problem. 

Gravitational clustering in an infinite space — static 
or expanding - starting from quasi-uniform initial con- 
ditions, is intrinsically a problem out of equilibrium. By 
"quasi-uniform" initial conditions we mean that the ini- 
tial state is a particle distribution — specified, we will as- 
sume, by a stochastic point process |l8| — which has rel- 
ative fluctuations at all scales, of small amplitude above 
the scale characteristic of the particle "granularity" and 
decaying at infinitely large scales ^. One of the most basic 



results (see e.g. |T^, and also the appendices to this 
paper) about self-gravitating systems, treated in a fluid 
limit, is that the amplitude of small fluctuations grows 
monotonically in time, in a way which is independant 
of the scale. This linearised treatment breaks down at 
any given scale when the relative fluctuation at the same 
scale becomes of order unity, signalling the onset of the 
"non-linear" phase of gravitational collapse of the mass 
in regions of the corresponding size. In an infinite space, 
in which the initial fluctuations are non-zero and finite 
at all scales, the collapse of larger and larger scales will 
continue ad infinitum. The system can therefore never 
reach a time independent state, and in particular it will 
never reach a thermodynamic equilibrium ^. One of the 
important results from numerical simulations of such sys- 
tems in the context of cosmology is, however, that the 
system nevertheless reaches a kind of scaling regime, in 
which the temporal evolution is equivalent to a rescal- 
ing of the spatial variables ^ . This spatio-temporal 
scaling relation is referred to as "self-similarity" ^. It is 
observed, however, only starting from a restricted class 
of simple initial conditions — we will describe these in 
further detail below — and in the specific Einstein de Sit- 
ter (EdS) expanding universe [ p^ . The range of initial 
conditions to which it applies has been a point of discus- 
sion in the literature, and theoretical explanations of it 
typically restrict it to quite a narrow range of such initial 
conditions, and strictly to the EdS expanding universe. 
To see whether this kind of simple behavior is reproduced 
in the system we study, is thus a first point of interest. 
It is in fact the primary focus of this paper. 

One comment needs to be made about the use of a 
static (Euclidean) space-time. The problem of bodies in- 
teracting by their mutual Newtonian self-gravity in the 
infinite volume limit, taken at constant mean density, is 
in fact ill defined: the force on a particle depends on how 
the limit is taken. In order to remove this ambiguity one 
adds a negative background to cancel the contribution of 
the mean density — the so-called "Jeans Swindle" (see 
e.g. 24 1). As discussed in ||2^, this is equivalent to tak- 
ing the limit symmetrically about each particle on which 
we calculate the total gravitational force^. Then only the 



^ We note that in ^ a treatment of infinite self-gravitating 
systems in the framework of equihbrium statistical mechanics 
is developed by considering a "dilute" infinite volume limit, in 
which Af — > oo and y — > oo at N/V^/^ = constant, where A'' 
is the number of particles and V is the volume (see also ]lo[ for 
a more recent discussion) . This is not the physically relevant 
limit for the problem treated here, as we consider the infinite 
volume limit taken at constant density, i.e., N ^ oo and — > oo 
at N/V = constant. In this discussed further below, 

the system is intrinsically time dependent and never reaches a 
thermodynamic equilibrium. 

3 It is also implicit in the phrase "quasi-uniform initial conditions 
in infinite space" that, as noted above, the infinite volume limit 
here is taken at constant particle density, rather than in the 
"dilute" limit studied in BlM- 



^ This does, of course, not mean that the instruments of equi- 
Ihriii rn s tatistical mechanics are completely irrelevant. Saslaw 
(see |l9| and references therein) notably has developed a treat- 
ment of gravitational clustering in an expanding universe which 
approximates it as a "quasi-equilbrium" in which the thermody- 
namic variables evolve adiabatically with the expansion of the 
Universe. Another more formal exploration of the usefulness of 
some standard equilibrium techniques can be found in |2o| . 

^ Note that this term is here used in a different sense to that com- 
monly ascribed to it in condensed matter physics. In this con- 
text "self-similarity" usually implies that the spatial correlations 
themselves have invariance properties under rescaling (see, e.g., 
tes]). This is not necessarily the case in the present context. 

^ See pf| for a very clear discussion of this issue. It is also shown 
here that addition of the negative backgound is equivalent to 
regularizing the problem with a cosmological constant. 
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fluctuations of the density field generate the gravitational 
force. In the context of cosmological expanding universe 
solutions, this "swindle" is unnecessary as the expansion 
absorbs the effect of the mean density, and the pertur- 
bations to the comoving particle trajectories are indeed 
sourced only by the fluctuations (see, e.g., Q). This 
modification does not necessarily make the gravitational 
force well defined in general: whether it is well defined 
depends on the nature of the fluctuations in the density 
field at large scales. For the case of the shuffled lattice 
(SL) considered here, we have studied in detail the prop- 
erties of the gravitational force in and shown the 
force to be well defined in the presence of the canceling 
background. 



limit^. Our initial conditions are similar, but not iden- 
tical, to those used in cosmological simulations of the 
formation of structure in the Universe. In this context 
the initial conditions are usually given by simple cubic 
lattices, perturbed by correlated displacements, with rel- 
ative displacements between nearest neighbor particles 
which are small . The displacements are generated in 
reciprocal space starting from an input power spectrum 
(PS), i.e., what is usually called the "structure factor" 
in condensed matter physics, specifying the desired the- 
oretical density fluctuations. 

In this paper we describe systematically basic results 
on gravitational dynamics starting from SL initial condi- 
tions. Our principal results are the following: 



Previous works in the same spirit as this [g7| ^ |2g[ | 
have treated primarily the very simplest initial condition 
one can envisage: Poisson distributed particles with no 
initial velocity. One of the basic results which has been 
emphasized in these works is the role of nearest neigh- 
bor interactions at early times in forming structures (see 
also [0), giving rise to non-linear density-density cor- 
relations which are then observed to be reproduced at 
larger and larger scales as time evolves. At the same 
time the effects of amplification at larger scales — de- 
scribed by the fiuid limit in which the granular structure 
of the matter is irrelevant — is observed. When trying to 
address the basic issue of the relative importance of these 
mechanisms, one runs into the limits imposed by the sim- 
ple initial conditions: in a Poisson distribution a single 
parameter — the particle density, or equivalently mean 
inter-particle distance — controls both the amplitude of 
fiuctuations and the "granularity" of the mass distribu- 
tion. This limitation is one of the major motivations for 
the different class of initial conditions we study in this 
work, developing further some initial analysis of this case 
in Ig^: we consider lattices subjected to small random 
displacements. In this case there are now two parame- 
ters, the inter-particle distance £ and the amplitude A 
of the "shuffling" . Given the scale free nature of gravity 
it is in fact only the dimensionless combination 6 = A/£ 
which is physically relevant (while in the case of Poisson 
initial conditions there is effectively no free adjustable 
parameter) . When the dynamics of the SL is treated in 
the fluid limit, as we will see, conflgurations with dif- 
ferent S may also be trivially related. In particular we 
can consider systems with different d which have differ- 
ent discreteness properties which are equivalent in terms 
of their fluid description. This allows us to understand 
notably the aspects of the evolution of the system which 
can be accounted for in a description of the dynamics in 
a fluid limit, and those which require the discreteness of 
the system to be explicitly taken into account. This is an 
important point as almost all existing analytic results on 
infinite self-gravitating systems are derived in this former 



• Evolution from these initial conditions converges, 
after a sufficient time, to a "self-similar" behavior, 
in which the two-point correlation function obeys a 
simple spatio-temporal scaling relation. The time 
dependence of the scaling ( i.e. the quantity anal- 
ogous to the dynamical exponent in out of equi- 
librium statistical mechanics) is in good agreement 
with that inferred from the linearized fluid approx- 
imation. This result is a generalization of what has 
been observed, for "redder" initial PS {P{k) ~ fc" 
with n < 1), in simulations in an EdS universe 




• Between the time at which the first non-linear cor- 
relations emerge in a given SL and the convergence 
to this "self-similar" behavior, there is a transient 
period of significant duration. During this time, 
the two-point correlation function already approx- 
imates well, at the observed non-linear scales, a 
spatio-temporal scaling relation, but in which the 
temporal evolution is faster than the asymptotic 
evolution. This behavior can be understood as an 
effect of discreteness, which leads to an initial "lag" 
of the temporal evolution at small scales. 

• Simulations with different particle numbers, but 
the same large scale fiuctuations (as characterized 
by the PS at small k), converge after a sufficient 
time, not only to the same functional form of the 
correlation function (with the self-similar behav- 
ior), but to the same amplitude. This is further 
evidence that it is indeed the common large scale 
fluctuations alone which determine the amplitudes 
of the correlations, which are thus independent of 
the discreteness scale £. At early times, however, we 



It is also a question which is very relevant in the context of cos- 
mology, as it concerns the understanding of the discreteness ef- 
fects in simulations of dark matter, which intrinsically limit their 
precision. These simulations treat the gravitational clustering of 
point "macro-particles" , which typically correspond to the order 
of 10^'^ dark matter particles. 
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see manifest difference between the systems, typi- 
cally again characterized as a "lag" of simulations 
with larger I (and smaller S) . 

• The non-linear correlations when they first develop 
are very well accounted for solely in terms of two- 
body correlations. This is naturally explained in 
terms of the central role of nearest neighbor in- 
teraction in the build-up of these first non-linear 
correlations. 

• This two-body phase extends to the time of onset 
of the spatio-temporal scaling, and thus the asymp- 
totic form of the correlation function is already es- 
tablished to a good approximation at this time. We 
briefly discuss the significance of this quite surpris- 
ing finding. 

The paper is organized as follows. In the next section 
we briefly define a SL distribution and introduce the main 
statistical quantities we use in the analysis and their esti- 
mators. We discuss the numerical simulations and their 



analyses in Sec. III. Finally in Sec. IV we summarize our 
mains results and conclusions, and briefly discuss some 
of the many open problems which remains for future in- 
vestigation. 



II. SHUFFLED LATTICES AND STATISTICAL 
QUANTITIES 

We firstly describe (Sec. J lI A| ) the class of initial condi- 
tions we study. In Sec jil B| we define the statistical quan- 
titi es we will use to characterize the correlations, and in 
Sec. II C we specify how we estimate these quantities in 
our simulations. 



A. Definition of a Shuffled Lattice 

We use the term SL to refer to the infinite point distri- 
bution obtained by randomly perturbing a perfect cubic 
lattice: each particle on the lattice, of lattice spacing £, is 
moved randomly ("shuffled") about its lattice site, each 
particle independently of all the others. A particle ini- 
tially at the lattice site R is thus at x(R) = R + u(R), 
where the random vectors u(R) are specified by the fac- 
torised joint probability density function 



7'[{u(R)}]=np(u(R)). 



(1) 



R 



The distribution is thus entirely specified by p(u), the 
probability density function (PDF) for the displacement 
of a single particle. 

In this paper we will study evolution from SL with the 



FIG. 1: Projection on the z — Q plane of a SL with 32 
particles and 5 — 0.177. Due to the random shuffling with the 
given PDF each lattice "chain" parallel to the z-axis projects 
onto a small square. 



following specific PDF 
p{u) = 




if ue hA,A]3 
otherwise. 



(2) 



Each particle is therefore moved randomly in a cube 
of side 2A centered on the corresponding lattice site 
(Fig. g). Taking A -> , at fixed £, one thus obtains 
a perfect lattice, while taking A ^ cx) at fixed £, one ob- 
tains an uncorrelated Poisson particle configuration [ p3| . 
Given Eq. (H), the shuffling parameter A^ also gives the 
variance of the shuffling, i.e.. 



A^ = y" d^uv?p{u). 



(3) 



Our SL configurations are therefore specified by two pa- 
rameters: the lattice constant € and the shuffling pa- 
rameter A. An alternative convenient characterization is 
given by £ and the adimensional ratio 5 = t^jt. We will 
refer to the latter as the normalized shuffling parameter. 
It is thus the square root of the variance of the shuffling 
in units of the lattice spacing. 

In what follows we consider not an infinite SL, but a 
finite SL of N particles in a cubic box of size L — N^^^i 
(see Fig. |l|). We will consider the specific case of a simple 
cubic lattice^ in which thus the mean number density of 



particles is thus N/L^ 



no 



Further we will assign 



to all particles the same mass m, so that the average mass 
density is simply po = totiq. 

In Table | we list the various relevant parameters of 
the SL considered as initial conditions of the N body 



* We will discuss in the conclusions section the importanc of this 
specific choice for the PDF. 
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Name 


JVV3 


L 


I 


A 


5 


m/m64 


SL64 


64 


1 


0.015625 


0.015625 


1 


1 


SL32 


32 


1 


0.03125 


0.0553 


0.177 


8 


SL24 


24 


1 


0.041667 


0.00359 


0.0861 


18.96 


SL16 


16 


1 


0.0625 


0.00195 


0.03125 


64 


SL128 


128 


2 


0.015625 


0.015625 


1 


1 



TABLE I: Details of the SL used as initial conditions in the 
simulations reported in this paper. A'' is the number of par- 
ticles, L is the box size, I the lattice constant and A (5) the 
(normalized) shuffling parameter. The mass m of the particles 
is expressed in unit of that in SL64, i.e., m64. In the units cho- 
sen, the mass density in all these systems is po = Nm/L^ — 1. 
Note that SL64 and SL128 are "more shuffled" than all the 
others (i.e. larger shuffling parameter) while SL16 is the one 
which is the closest to a perfect cubic lattice. Note that SL128 
differs only from SL64 by the size of the box. 

simulations (NBS) which we report here. We will explain 
below the criteria used for these choices. 



B. Statistical Characterization of Correlation 
Properties 

The microscopic number density function for any par- 
ticle distribution is given by 

TV 

n(x) =^(5d(x-x,) , (4) 

1=1 

where is the position of the z-th particle, Sd is the 
Dirac delta function and the sum is over the N particles 
of the system. 



It is useful also to note that one can write 

= (8) 
m 

where (n(r))p, the conditional average density, is the (en- 
semble) average density of points in an infinitesimal shell 
at distance r from an occupied point ^ . We will make use 
of this relation is estimating ^(r) below. 

In the evolved self-gravitating systems we study ^(r) 
will invariably be a monotonically decreasing function of 
r. It is then natural to define a scale A by 

C(A) = 1 (9) 

which separate the regime of weak correlations (i.e. 
^(r) ^ 1) from the regime of strong correlations (i.e. 
C(r) 3> 1). In the context of gravity these are what are 
referred to as the linear and non-linear regimes, as a lin- 
earized treatment of the evolution of density fiuctuations 
is expected to be valid in the former case. Given the 
form of Eq. (^) it is clear that A as defined by Eq. (^) is 
an appropriate definition of the homogeneity scale of the 
system. This scale gives then the typical size of strongly 
clustered regions. 

The exact analytic expression for ^(r) in a SL is given 
in . We do not reproduce it here as it is a complicated 
expression, which we will not in fact make use of. In our 
case, as in the case of a perfect lattice and a Poisson dis- 
tribution (which, as we have noted correspond to specific 
limits of the SL) ^(r) < 1 everywhere: there is no strong 
clustering. In such a situation the homogeneity scale is of 
order of the average distance between nearest neighbors 
(NN), which we will denote by A. Thus when we refer to 
the homogeneity scale we will mean A in absence of non- 
linear clustering and the scale given by Eq.(p[) otherwise. 



The two-point correlation function 



2. The mass variance 



For a system such as we consider here, in which the 
mean density is well defined and non-zero, it is convenient 
to define the density contrast: 



<5(x) = 



7l(x) - HQ 



HQ 



(5) 



In order to characterize the two-point correlation prop- 
erties of the density fluctuations, one can then use the 
reduced two-point correlation function: 



|(r) = (<5(x + r)5(x)) 



(6) 



where (...) is an ensemble average, i.e., an average over all 
possible realizations of the system. In a distribution of 
discrete particles ^(r) always has a Dirac delta function 
singularity at r = 0, which it is convenient to separate 
by defining ^(r) for r 7^ (the "off-diagonal" part) [Ssl: 



|(r) = -5^(r)+e(r) . 
no 



(7) 



For particle distributions with a well defined average 
density it useful also to consider an integrated quantity 
such as the normalized variance of particle number (or 
mass) in spheres, defined as : 



a\r) 



{N^r)) - {N{r)y^ 
WW' 



(10) 



where N{r) is the number of particles in a sphere of ra- 
dius r. Then cr'^{r) can be used, in a manner similar to 
that described above for ^(r), to distinguish a regime of 
large fluctuations from a regime of small fluctuations. It 



^ For the more general case of non-uniform distributions, such as 
fractal particle distributions [ p3[ , in which no is zero, this is the 
basic statistical quantity for the characterization of two-point 
correlation properties, rather than ^(r) (which is then not de- 
fined) . 
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is simple to find the explicit expression for a'^(r), which 
gives it as a double integral of ^(r) Q. 

One can show that the normalized variance in real- 
space spheres defined in Eq. behaves in a SL as 
CT^(r) cx r~'* at large r, compared to cr^(r) oc in a 



Poisson distribution 



The behavior cr {r) cx 



is in fact the fastest possible decay of this quantity [g3| . 
This means that the SL belongs to the class of distribu- 
tions which may be termed super-homogeneous [|4| (or 
hyper-uniform [^). Such systems have mass fluctua- 
tions which are depressed with respect to those in an 
uncorrelated Poisson distribution. 



3. The Power Spectrum 

Since we consider distributions which are periodic in 
a cube of side L, we can write the density contrast as a 
Fourier series: 



(5(x) ^ exp(ik • x) (5(k) 



(11) 



with k e {(27r/L)n| n e I?]. The coefficients 5(k) are 
given by 

(5(k) = / (5(x) exp(-ik • x) d^x . (12) 

The PS of a particle distribution^" is then defined as (see 
e.g. MM) 



F(k) = -(|5»P) 



(13) 



In distributions which are statistically homogeneous, 
which is the case here^-'^, the PS and reduced two-point 
correlation function ^(r) are a Fourier conjugate pair. 

The exact expression for the PS of a SL is simple to 
derive. One finds (see ^ ) 



P(k) = 



l-|p(k)P 



L3^|p(k)p5K(k,n^) (14) 



where p(k) is the Fourier transform of the PDF for the 
displacements p(u) (i.e. its characteristic function), and 
i5k is the three-dimensional Kronecker symbol. For the 
specific p(u) given in Eq. (||) we have 



\m?= n 



(fc,A)2 



(15) 



Inserting this expression in Eq. (|lj) one obtains an exact 
explicit analytic expression for the PS of a SL in terms of 
the two parameters £ and A. It is simple to verify that 
taking A/£ = 6 ^ oo, at fixed £, one obtains P{k) = l/no 
(as expected, since one obtains in this limit a Poisson 
distribution). Further one always approaches this same 
behavior (as required [^Sj) in the limit k —^ oo. Further, 
at small k (i.e. fc <C 27r/£), we obtain 



P(k) 



IkpA^ 
3no 



(16) 



We note that this result can actually be found (see |3^, 
|3^ ) directly from Eq. (p^), without assuming a specific 
form for p{u). One need only assume that p{u) has a 
finite variance, equal to A^. Thus the small k behavior 
of the PS of the SL does not depend on the details of the 
chosen PDF for the displacements, but only on its (finite) 

1 9 

variance. . 

Finally note that the mass variance can actually be ex- 
pressed simply as an integral in reciprocal space of the PS 
multiplied by an appropriately normalized Fourier trans- 
form W{k]r) of the spherical window function, being 
outside the sphere and 1 inside it |33]: 



(r) 



1 



(27r)3 



d^kP{k)\W{k;r)\^' 



(17) 



4- The nearest neighbor distribution 

A very useful and simple statistical quantity which 
characterizes small-scale clustering properties of a parti- 
cle distribution is the nearest neighbor (NN) PDF uj{r). 
It is the probability density for the distance between a 
particle and its NN Q, i.e., io(r)dr gives the probability 
that a particle has its NN at a distance in [r, r -\- dr]. If 
one neglects correlations of any order higher than two, it 
is simple to show that w(r) is related to the conditional 
density {n{r))p (and thus to ^(r), given Eq. (@)) through 

13 



w(r) dr = 1^1 - y Lu{s) dsj ■ 47rr^(n(r))pdr , (18) 

This relation will be very useful to us here because it is 
valid in particular when clustering is dominated, at early 
times, by individual pairs of particles falling toward each 
other. 



We use here the term for this quantity commonly employed in 
cosmology, rather than "structure factor" which is more habitual 
in the context of condensed matter and statistical physics. Note 
also the normalisation, which corresponds to P(k oo) — » 
rather than unity. 

For the lattice and SL the ensemble average is defined over the 
set of lattices rigidly translated by an arbitrary vector in the unit 
cell. 



Note that the behavior limk^o P(k) = is an equivalent way 
of stating the property of super-homogeneity of the distribution 

relation follows if one assumes that the probability of finding 
a particle in [r, r -\- dr] given that there is a point at r = is the 
same whether the condition that there be no other point in [0, r] 
is imposed or not. 
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C. Estimation of Statistical Quantities 

In order to estimate P(k) and ^(r) in a given particle 
configuration, i.e., in a single realization of the evolved 
SL, we calculate averages in spherical shells in real or 
reciprocal space. This means that we consider only the 
dependence of these quantities on the modulus of their 
arguments and we will therefore use the notation P{k) 
and ^(r) in the rest of the paper. 

The PS is obtained from (5(k) by means of the rela- 



tion 
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Pik) 



N(k) , ^ 

^ ' fc<|k'|<fe+(SA; 



|^(k')|^ 



(19) 



where iV(fc) is simply the number of vectors k' considered 
in the sum. Note that to speed up the calculations, not all 
the vectors k' for a given modulus are taken into account: 
at large k the density of vectors considered is smaller than 
at small k. 

The function ^(r) is estimated by first calculating 
(n(r))p (see Eq. ||) Q. As already mentioned the lat- 
ter gives the average density in a spherical shell of radius 
r, and thickness <C r, centered on an occupied point. 
Thus we estimate it as 



1 



E 



(20) 



where Ni{r) is the number of particles in the spherical 
shell of radii r, r + 6r, volume V{r, Sr), centered on the 
i*^ particle of a subset of A^c < particles randomly 
chosen among the N particles of the system. 
The mass variance can be simply estimated by 



1 



1=1 



{Nir))f (21) 



where Ni(r) is the number of particles contained in the 
i*^ (with i — l..Nc) randomly placed sphere of radius r 
and {N{r)) its average. 

Finally the NN distribution a;(r) is computed directly 
by pair counting. 



which is based on a tree algorithm for the calculation of 
the force, allows one to perform simulations in an infinite 
space, using the Ewald summation method The po- 
tential used is exactly equal to the Newtonian potential 
for separations greater than the softening length e, and 
regularized at smaller scales. For what concerns the in- 
tegration parameter we have performed several tests to 
check the stability of the results at the level of numerical 
precision we consider in this work 

We have considered as initial conditions the set of five 
SL described in Table ||. We now explain the reasons for 
our choices of the parameters given. 

Firstly it is important to note that, in the limit of the 
pure (i.e. un-softened) gravitational evolution of an in- 
finite SL, there is only one parameter which can change 
the dynamical evolution non-trivially . This is S, the 
normalized shuffling parameter (i.e. normalized in units 
of the lattice spacing £) . Because gravity has no preferred 
length scale the gravitational dynamics of two infinite SL 
with the same S, but different lattice spacing £, can be 
trivially related: a rescaling of length scales is equivalent 
to a rescaling of time, so that the configurations of one 
can be mapped at any time onto the configuration of the 
other at a different time. The same is true for changes 
of the mass of the particles: two SL with the same S, 
the same £, but different particle masses, are related by 
a simple scaling of the time variable. Indeed any two SL 
with the same 5, are strictly equivalent to one another in 
time if they are related to one another by any simultane- 
ous scaling of £ and the particle mass which leaves their 
mass density po fixed. 

For softened gravity in a finite box, the same length 
scale transformations can relate trivially different SL 
with the same d. In this case the relevant parameters 
to distinguish two SL evolved in these simulations are 
thus three, which we can take as S, £/L = N^^"^ and e/L. 

We have chosen our (arbitrary) units of length, mass 
and time as follows. Our unit of length is given by the 
box side of the SL64 simulation and our unit of mass 
by the particle mass in this same simulation. A natural 
choice for the unit of time is the so called dynamical time, 
defined as 



III. GRAVITATIONAL CLUSTERING IN A 
SHUFFLED LATTICE: RESULTS FROM 
NUMERICAL SIMULATIONS 



1 



Tdyn 



%/47rGpo 



(22) 



As unit of time here we have made a slightly different 



Details of Numerical Simulations 



We have performed a set of numerical simulations us- In order to test the numerical accuracy of the simulations we have 

ing the freely available code Gadget p^, p3 . This code, also compared the early times evolution with the prediction of 

the linearised treatment of the early time evolution, as described 
in detail in jio) . 

By "dynamical evolution" we mean the ensemble average prop- 
erties of the clustering etc.. We thus suppose that this average 
For simplicity in this paper we use the same symbol for the en- is recovered in a single realization of an infinite volume system, 

semble average quantity and for its estimator. i.e., that spatial ergodicity applies. The unaveraged dynamical 

We use periodic boundary conditions in this estimation (as in evolution will of course vary in detail from one realization to 

the simulations). another. 
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choice of the pre-factor, with Tdyn = 1.092 

In the "reference" SL64 simulation we have chosen 
6 = 1. Our choices of the parameters for the other five 
simulations can be understood as follows: 

• The particle masses are chosen so that the mass 
density is constant. Thus the dynamical time is 
the same in all simulations, which is convenient for 
comparison, as we will see, as this is the unique 
time scale of these systems in the fluid limit. 

• The softening e is the same in all the simulations. 
We have chosen e = 0.00175 in our units, which 
means that it is, in all the simulations, significantly 
smaller (at least a factor of ten) than Aq, the initial 
average distance between NN 

• The box size is the same in all but one simulation. 
This latter simulation (SL128), which is the biggest 
one, is used to test the accuracy with which our re- 
sults are representative of the infinite volume limit 
(at fixed mass and particle density). Thus it is cho- 
sen to have the same parameters to SL64, differing 
only in its volume (which increases by a factor of 



For each of the four other SL simulations we change 
the number of particles N, which fixes £. We have 
then chosen S so that the PS at small k has the 
same amplitude. From Eq. ( p^ ) it it easy to see 
that this requires, in our chosen units. 



no 



SL64 



64" 



(23) 



The PS of the SL described in Table | are shown in 
Fig. ||. We see, up to statistical fluctuations, that the 
spectra are indeed of the same amplitude at small k. 
Note that the Nyquist frequency fcjv = tt/^ in fc-space 
translates to the right with increasing particle number. 

The particles are assigned zero velocity in the initial 
conditions (at t = 0), and, as has been underlined, the 
simulations are performed in a static Euclidean universe, 
i.e., without expansion or non-trivial spatial curvature. 
We have run the simulations SL16, SL24, SL32 and SL64 
up to about time 6 and the SL128 up to time 8 as for 
longer times the simulations begin to be dominated by a 
single non-linear structure, a regime in which we are not 
interested since it is evidently strongly affected by finite 
size effects. 



This corresponds to time in units of 1000 seconds for a mass 
density of Ig.cm"^. 

The smallest value of Aq is that in SL64 where it is equal to 
0.55/64 0.0086 as in a Poisson distribution with the same 
number density |p^. 
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FIG. 2: The PS (averaged in spherical shells) of the SL con- 
figurations specified above in Tab. | as a function of the mod- 
ulus of k. The solid line is the theoretical (oc k^) behavior 
for small k given by Eq. (p^. At large k, the four PS are 
equal to 1/no, with the correspond ing value of no. The peaks 
arise from the second term in Eq. (|l4[). The four arrows show 
the different Nyquist wave-numbers multiplied by two for the 
SL configurations in order of increasing number density from 
left to right: this corresponds to the expected location of the 
first peak in each case. Note that, as we have discussed in 
Sect. [IC, not all the vectors k are considered in the estima- 



tion of the PS and therefore only a subset of all the peaks is 
detected (each peak corresponds to a very narrow band of k 
so it can be easily missed). 



B. Results 

In this section we analyze the results of the numeri- 
cal evolution of the SL described in Table | in terms of 
the statistical quantities discussed above. In the first two 
subsections we restrict ourselves to the study of the evo- 
lution of SL128, i.e., the largest simulation we have run. 
This is then our reference point with which we compare 
the evolution of the other initial conditions. 



1. Evolution of the Power Spectrum 

The evolution of the PS in SL128 estimated by using 
Eq. ( p9| ) is shown in Fig. ^ Along with the numerical 
results is shown the prediction for the evolution of the 
PS given by the linearized fluid theory (see App. ^: 



P(k,i) =P(k,0)cosh2 (Vrdyn) 



(24) 



We observe that: 



The linear theory prediction describes the evolu- 
tion very accurately in a range k < k*{t), where 
k*{t) is a wave- number which decreases as a func- 
tion of time. This is precisely the qualitative be- 
havior expected as linear theory is expected to hold 
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FIG. 3: Evolution of the PS in SL128 (solid lines — label 
FG): the curves are for time equal to 0,2,4,6,8 (from bottom 
to up). The dashed lines labeled with LT show the predictions 
of fluid linear theory, i.e., Eq. ( |24| ) with P(k,0) measured in 
the simulation at t = for the same time steps. The arrow 
labeled "fcjv" shows the value of the corresponding Nyquist 
frequency fcjv = n/i. 



only above a scale which, in real space, increases 
with time, and, in reciprocal space, decreases with 
time. We note that at t — 6 only the very smallest 
/c-modes in the box arc still in this linear regime, 
while at i = 8 this is no longer true. We will discuss 
below a more precise quantification of the validity 
of the linearized approximation. 

• At very large wave-numbers (fc > fcjv) the PS re- 
mains equal to its initial value 1 /no . This is simply 
a reflection of the necessary presence of shot noise 
fluctuations at small scales due to the particle na- 
ture of the distribution. We note that the value 
of k at which this behavior is attained increases 
somewhat from its initial value and then remains 
roughly stable. We will comment further on the 
significance of this fact below. 

• In the intermediate range of A:, i.e., k*{t) < fc < k^, 
the evolution is quite different, and slower, than 
that given by linear theory. This is the regime of 
non-linear clustering. 

These results concerning the validity of linear theory 
at sufficiently small fc, and in a range which decreases 
as a function of time, are completely in line with what 
is observed in cosmological simulations, in an expand- 
ing universe (see e.g. pd|). In this context simulations 
typically start from lattices with correlated perturba- 
tions representing spectra which are much "redder" than 
P{k) ~ k^, typically P{k) - /c" with -3 < n < 0. That 
the same behavior is seen in a static universe for this 
"bluer" PS is, however, expected. Indeed, on the basis 
of simple considerations (see e.g. W^) — which do not 



FIG. 4: Behavior of the absolute value of the correlation func- 
tion |C(r)| in SL128 at times t = 0, 1, 2, 4, 6, 8. Note that val- 
ues such that |^(r)| < 0.01 — 0.1 are below the level of noise of 
the estimator estimated by using the normalized variance in 
spherical shells (see text) : This gives a limit below which the 
noise in the estimator dominates over the signal. The arrows 
shows the value of the lattice spacing £ and the initial aver- 
age distance between nearest particles Ao, while the dotted 
vertical line corresponds to the smoothing e. 



make use of the expansion of the universe — about the 
long-wavelength (i.e. small k) perturbations generated 
by non-linear motions on small scales, one anticipates 
that linear theory should be valid at small k for any ini- 
tial PS with P{k) ^ fc" and n < 4. The reason is that 
such non-linear motions, which preserve locally mass and 
center of mass, can generate at most a PS at small k with 
the behavior P{k) ~ k^. 



2. Evolution of the Two-Point Correlation Function 

We consider now the evolution of clustering in real 
space, as characterized by the reduced correlation func- 
tion ^(r). We focus again on SL128. In Fig. ^ is shown 
the evolution of the absolute value \^{r)\ in a log- log plot. 
In the figure is shown also, for comparison, at large scales, 
the level of the typical fiuctuations expected in the esti- 
mator of ^(r) This indicates that, at larger separa- 
tions, the noise in the estimator is expected to dominate 
over any underlying physical correlation which may be 
present. 

We observe that: 



This estimate is obtained by assuming that the variance in the 
shells employed in the estimator decay as ^ ""shell ('')/-^<^ oc 
where a'^jjell('') variajjce in shells, defined analogously to 

that in spheres [cf. Eq. (|lO|)] , and Nc is the number of centers 
used to calculate ^(r) [cf. Eq. E^]. 



10 



< 0.006 




012345678 



FIG. 5: Evolution in time of A{t), the average distance be- 
tween nearest neighbors, in SL128. It decreases at early times 
and then stabilizes at A « 2e. 

• Starting from ^(r) ^ 1 everywhere, non-linear cor- 
relations (i.e. ^(r) ^ 1 ) develop first at scales 
smaller than the initial inter-particle distance Aq. 

• After two dynamical times the clustering develops 
little at scales below e. The clustering at these 
scales is characterized by an approximate "plateau" 
at £^{r) ~ 10^. This stabilization of the system at 
small scales is also evident in Fig. ||, which shows 
the evolution of the mean distance between nearest 
neighbor particles as a function of time. The stabi- 
lization in time of the scale in k space at which the 
PS reaches its asymptotic (constant) value, which 
we observed above, is just the manifestation in re- 
ciprocal space of this same behavior. 

• At scales larger than e the correlations grow con- 
tinuously in time at all scales, with the scale of 
non-linearity [which can be defined, as discussed 
above, by ^(A) = 1] moving to larger scales. 

From Fig.^ it appears that, once significant non-linear 
correlations are formed, the evolution of the correlation 
function ^(r) can be described, approximately, by a sim- 
ple "translation" in time. This suggests that £,{r,t) may 
satisfy in this regime a spatio-temporal scaling relation: 

ar,t)^E{r/R,{t)) , (25) 

where Rs{t) is a time dependent length scale which we 
discuss in what follows. In order to see how well such 
an ansatz describes the evolution, we show in Fig. ^ an 
appropriate "collapse plot": £,{r,t) at different times is 
represented with a rescaling of the x-axis by a (time- 
dependent) factor chosen to superimpose it as closely as 
possible over itself at t = 1, which is the time from which 
the "translation" appears to first become a good approx- 
imation. We can conclude clearly from Fig. ^ that the 
relation (p5|) indeed describes very well the evolution, 




r 



FIG. 6: Collapse plot of £,{r,t): for each time t > 1 we have 
rescaled the a;-axis by a time-dependent factor to collapse all 
the curves (dashed ones) to that at time f = 1. We have added 
for comparison ^{r,t — 8) without rescaling ("w. resc", con- 
tinuous line). 
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FIG. 7: Evolution of the function Rs{t) in SL128 (points) 
compared with its prediction ("expected" and "theoretical") 
for different values of the time scale of onset of self-similarity 
(see text for details). Both lines corresponds to Rs{t) oc 
exp[(2/5)i/Tdyn]- Also shown is the corresponding prediction 
for Poisson initial conditions, -Rs(t) oc exp[(2/3)t/rdyn]- 



down to separations of order e, and up to scales at which 
the noise dominates the estimator. 

In Fig. is shown the evolution of the rescaling factor 
Rs{t) found in constructing Fig. o, as a function of time, 
with the (arbitrary) choice Rs{l) — 1. Shown in Fig. 
are also three (theoretical) curves, which we will explain 
in the next subsection. Before this we remark on two 
further aspects of the scaling relation which are worth 
noting: 

• The function S(r), when it is larger than ~ 0.1, can 
be well approximated by a simple power law with 
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time), is referred to as self-similarity. Such behavior is 
expected in an evolving self-gravitating system (see e.g. 
]l6| , |3^ ) because of the scale- free nature of gravity, if 
the expanding universe model and the initial conditions 
contain no characteristic scales. Initial conditions in N- 
body simulations do, however, necessarily contain one 
such scale, which is associated to the particle discrete- 
ness (i.e. the grid spacing £ in the case of a perturbed 
lattice). Further, as we have discussed above, simulations 
introduce (at least) two further scales: the box-side L and 
force softening e. Thus self-similarity is expected to be 
observed in N body simulations of an Einstein-de Sitter 
model (i.e. a flat matter dominated universe), starting 
from pure power-law initial PS P{k) ^ fc", if all effects 
associated with these length scales can be neglected. 



FIG. 8: Comparison of (,{r,t) measured in SL128 with the 
formula in Eq. ([26[), using the rescaling of Eq. (^^. For clarity 
the amplitudes of the different curves have been rescaled (by 
a factor 3'~^). Moreover for all the curves, we have plotted 
only the scales such that £,{r,t) ^0.1 since Eq. (Eq) does 
not describe smaller amplitude correlations. The vertical line 
corresponds to the softening length e. 



an exponential cut-off: 



S(r) « A 



R 



exp 



(26) 



where we have estimated (see Fig. ^) the following 
values for the three parameters: A = 40, 7 = 0.28 
and R = 1.45 • 10~^. This last parameter gives the 
normalization of £,{r,t) at t — 1. In order to see 
how well this fit describes the evolution, we show 
in Fig. 1^ both the data and the curves inferred from 
it, using Rs{t) as measured. 

• Since we have defined the homogeneity scale A by 
^(A) = 1 it is clear that, once the spatio-temporal 
scaling relation is valid, we have X{t) cx Rs{t). 

• Since the PS and mass variance are simply related 
to ^(r), we expect the scaling relation to be re- 
flected in one for these quantities as well. We will 
see to what extent this is the case below. 



3. Spatio-temporal scaling and "self-similarity" 

We have observed that, from t ^ Tdyn, the two-point 
correlation function, at least down to ^(r) ^ 0.01 (level 
of estimator noise) obeys to a good approximation the 
spatio-temporal scaling relation Eq. (P5[), with the mea- 
sured Rs{t) shown in Fig. ^ In this section we discuss 
this result, in particular its relation to similar behaviors 
which have been studied in cosmology. 

In the context of cosmological N body simulations this 
kind of behavior, when Rs{t) is itself a power law (in 



On theoretical grounds there are different expectations 
([pH |3^ , Q) about the range of exponents n of the PS 
which should give self-similar behavior. Effects coming 
from the particle discreteness are expected to become less 
important as the PS becomes "redder" (i.e. smaller n, 
with more relative power at larger scales), while a PS 
which is too "red" will become sensitive to finite size ef- 
fects (i.e. to the box size). A more quantitative analysis 
of the dependence of dynamically relevant quantities (e.g. 
the variance of velocity and force fields) on these ultra- 
violet and infra-red cut-offs suggests that self-similarity 
should apply in the range — 1 < n < -1-1, and such be- 
havior has in fact been observed, to a good approxima- 
tion, to apply in simulations in an EdS universe of such 
spectra ^1, |2^. While there has been considerable dis- 
cussion also of the case — 3<n<— lin the literature, 
with different conclusions about the observed degree of 
self-similarity, the case n > 1 has remained open In 
our discussion below we will see in greater detail why the 
cases n > 1 and n < 1 are expected to be possibly very 
different with respect to "self-similarity" . 

Our results above clearly suggest that what we have 
observed is a simple generalization of this self-similarity 
to a static universe (in which there is evidently also no 
characteristic length scale), and to the case n = 2. Let us 
examine more carefully whether this is the case, by gen- 
eralizing to the static case the argument (see e.g. [ |l6| ) 
used to derive the power-law behavior of Rg (t) in an ex- 
panding universe. 

In order to derive this behavior of Rsit), we assume 
that the spatio-temporal scaling relation holds exactly, 
i.e., at all scales, from, say, a time ts > 0. For t > ts we 



The reason why the case n > 1 has not been studied numerically 
appears to be twofold: firstly, it is not of direct interest to "real" 
cosmological models which typically describe PS with exponents 
in the range —3 < n < —1; secondly, the simulation of such 
initial conditions is considered "hard to simulate" (see e.g. [p^). 
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have then 

P(k,t)=/ exp(-ik-r)^(r,i)d^r 

= Rl{t)( exp(-ii?,(i)k-x)S(|x|) d^x (2^) 
= Rl{t)P{R,{t)]^,U) . 

where we have chosen Rs{ts) = 1- Assuming now that 
the PS at small k is amplified as given by linear theory, 
i.e., as in Eq. (|^, one infers for any PS P{k) ~ fc" (and 
n < 4 so that linear theory applies): 

R m I r 2(^-t.) 1 

yCOSh;;^y L(3 + n)Tdvn 

(28) 

In the asymptotic behavior the relative rescaling in space 
for any two times becomes a function only of the differ- 
ence in time between them so that we can write 

e(M + Ai)=e(^,^) ; R.m-e^^^. 

(29) 

This is analogous to what is called self-similarity in EdS 
cosmology. In that case the linear theory describes a 
growing and a decaying mode, both of them power laws in 
time. Asymptotically Rs{t) is thus itself a simple power 
law 

Let us now return to Fig.^ In addition to the measured 
values of Rs (t) the figure shows two curves corresponding 
to Eq. (H) with n = 2. The first (labeled "expected") 
corresponds to taking = 1 in the derivation above, i.e., 
assuming that the scaling relation holds at all scales for 
t > 1. The second (labeled "expected (resc)") is the same 
functional behavior, but rescaled by a constant to give a 
good fit to the larger time (from t > tg ^ 2.5) behavior. 
This latter behavior is clearly very consistent with the 
relation given in Eq. (p8|): starting from this time the 
slope is very close to constant and equal to | in units of 

Our results are thus indeed clearly well interpreted as a 
generalization in a static universe of self-similarity as ob- 
served in simulations EdS universes, for "redder" spectra. 
This self-similarity sets in, however, from about ts — 2.5, 
while we observed the spatio-temporal scaling relation 
already to apply approximately from t ~ 1 (w ''dyn)- 
Note that the fact that the functional behavior of Rs{t) 
in 1 < t < 2.5 is inconsistent with Eq. (|2^ ) with n = 2 
implies that the spatio-temporal scaling relation cannot 
hold at all scales at these times: specifically it cannot 
hold at small fc, where P(k) (x fc^, as we have seen that 
at these scales the PS is linearly amplified at this time. 



One has |l6[ ^ P(k,t) oc t^/^fc" at small k, and thus Rsit) oc 

t3(3 + ") . 



A possible explanation for this behavior is suggested by 
the third curve (labeled "Poisson" ) shown in Fig.^. This 
curve corresponds to Eq. ( p^ ) with n = and tg — I. 
The fact that it fits the points reasonably well — al- 
though not so well as the n — 2 theoretical curve for 
t > 2.5 — suggests the following interpretation: be- 
tween 1 < t < 2.5 we are observing a first phase of 
"self-similarity", restricted to smaller scales, where the 
initial PS is roughly flat (i.e. Poisson-like with n — 0) 
in a small range of k around the Nyquist frequency (see 
Fig. ^ . Such an interpretation is consistent with the fact 
that the ^(r) in the non-linear regime observed in sim- 
ulations from Poissonian initial conditions is, to a very 
good approximation, the same as that observed from SL 
initial conditions |^ , On the other hand, the wave 
modes at which the PS is Poisson-like are very large — 
of the order of the inverse of the inter-particle spacing — 
and so the observation of apparent self-similarity driven 
by these fluctuations is somewhat surprising: such be- 
havior is expected, as we have discussed above, when 
the effects of discreteness may be neglected. We will 
see below that this interpretation of the spatio-temporal 
scaling observed in the correlation function at non-linear 
scales at early times — as a first phase of self-similarity 
driven by Poisson fluctuations at small scales — is not 
correct. In particular it is reproduced in the smaller sim- 
ulations we will analyze below in which there is no initial 
Poisson plateau around k^. Further we will see that 
the form of the non-linear correlation function is already 
established at times when two-body correlations due to 
nearest neighbor interactions are the dominant source of 
correlation at these scales. 



4- Evolution of the mass variance 

In this section we study the normalized mass variance 
0-2 (r), deflned in Eq. (pX|). Through the study of this 
quantity we can probe further the scaling properties (and 
self-similarity) we have just seen. We can also explain 
and see the interesting and non-trivial differences in this 
respect between the case of a PS with n < 1 and n > 1. 

Given that the mass variance is expressible (cf. Eq. p^ ) 
as an integral of the PS, one might anticipate that it will 
show, at large r, the same behavior as the PS at small 
k, i.e., we expect to flnd the simple scale independent 
amplification of linear theory: 

a2(r,t) = A,2(i)a2(r,0) (30) 

where 

A„2{t) = cosh2(t/Tdy„) cx i?^+"(i) (31) 

for a PS P{k) cx fc" at smaU k. 

In Fig. ^ is shown the temporal evolution in SL128 of 
cr2(r). At each time we observe at large r the behavior 
0-2 (r) cx 1 /r"* characteristic of a SL. The dotted lines show 
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FIG. 9: Evolution of the mass variance in SL128 at times t — 
0,2,4,6,8, together with the predictions of Eq.(p2[) (labeled 
as LT). 
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FIG. 10: Behavior of P(A;max, i)/fcmax as a function of time 
measured in SL128. The dashed line represents the behavior 
given in Eq. (ba). 



the best fit to the behavior of Eq. (|30|) above, which we 
find is 

A,2{t) = cosh«/^(t/Tdy„) oc . (32) 

rather than the anticipated behavior of Eq. (|l|) for n = 
2. 

The reason for this discrepancy is, as we now discuss, 
very simple. It is of importance as it is makes explicit 
the difference between the cases of PS with n < 1 and 
n > 1. Indeed examining the integral Eq. (p^ ) in closer 
detail it turns out that there is a qualitatively different 
behavior in the two cases. For —3 < n < +1 the integral 



is dominated by modes k • 



and one has 



a'^[r,t) w Ck^P{k,t) 



(33) 



where k = 1/r and C is a constant pre- factor which de- 
pends on n. From this it follows that linear amplification 



FIG. 11: Behavior of A(a, t), the length for which cr^(r, t) = a 
for a = 0.1, 1.0, 10.0, 20.0 in the simulations SL128. 



of the PS at small k gives linear amplification of the mass 
variance at large scales. For n > 1, however, the integral 
in Eq. ( p^ ) with P{k) ~ fc" at all k diverges, and an ul- 
traviolet cut-off kc above which P{k)/k decays to zero is 
required to regulate it . The effect of the cut-off is to 
give 

a^r) ^ k-^P{k,)/r^ , 

at sufficiently large r. Thus, for n > 1 the evolution of the 
mass variance at large r (and thus at small amplitude) 
is sensitive to the evolution of the cut-off in the PS (and 
the amplitude of the PS at this cut-off). From FigJ| we 
expect that in our system the role of fee will be played 
by A:max(i), the wavenumber at which the PS reaches its 
maximum, and so we will have, at large r 



a\r,t) 



max 

it)P{kme.^,t) 



(34) 



From FigJI we see that fcmax is clearly in the range in 
which the amplification in k space is non-linear. Thus 
the evolution of this quantity, even at very large scales, 
is determined by modes in k space which are in the non- 
linear regime. Given the time evolution we have observed 
for (r, t) in Fig.^, we must have 



(35) 



In Fig. |lO| we see that this behavior is indeed well ap- 
proximated. It is in fact evidently the direct consequence 
of the self-similarity as it is reflected in the variance, and, 
equivalently, in k space. To the extent that both quan- 
tities approximate the self-similarity observed in ^(r, i). 



Such a cut-off necessarily exists in any particle distribution as 
P{k) cannot diverge for large k. One necessarily has, as we have 



discussed, P{k) — > 1/no for A; — > oo [p3|. 
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any length scale derived from either the variance or PS 
must scale oc Rs{t). Thus, in particular, k^]^^ (x Rs(t) 
and the maximum value of the PS, which has dimen- 
sions of volume, must scale as P(fcniax,i) oc Rl{t) in this 
regime. 

It is instructive also to examine a little further how 
the spatio-temporal scaling behavior, and self-similarity, 
are approximated in the variance. In order to illustrate 
this we consider the temporal evolution of scales X{a, t) 
defined by the relation 



(36) 



where a is a chosen constant. If there is a spatio-temporal 
scaling in the system we should find that A(a, t) cx Rs(t). 
In particular any choice a ~ 1 gives, as we discussed in 
Sec.|l[ a reasonable definition for the homogeneity scale, 
which should be equivalent to the one we have taken 
above {£,{X,t) = 1 once non- linear clustering has devel- 
oped). In Fig. y we show A(a, t) for a = 20, 10, 1, 0.1; 
also shown are curves proportional to Rs{t) in the self- 
similar regime (i.e. as given by Eq. (^ ) with n = 2). 
The figure illustrates nicely how the scaling applies only 
at large scales (corresponding to smaller fluctuations) ini- 
tially and then propagates to smaller (more non-linear) 
scales. At the time t « 2.5, which we identified above 
in our analysis of ^(r, t) as the time from which self- 
similarity is well approximated, the scaling behavior 
given by Rs(t) is manifestly well approximated well for 
a ^ 1, i.e., into the non-linear regime. We do not, how- 
ever, see a behavior consistent with the hypothesis that 
the evolution prior to this time (1 < t < 2.5) is self- 
similar and associated to a PS with n = around kj^: 
this would correspond as in Fig.^ to a faster evolution of 
the scales shown here at these times, which is not what 
is observed. 



Self-similarity and the regime of validity of linear theory 



The derivation of Rs{t) in Eq. (^8|) explains implic- 
itly the physical origin of the self-similar behavior: if the 
small k PS is a simple power law, the evolution of the 
two-point correlation function is self-similar, with Rs{t) 
given by Eq. (p8|), in the approximation that fluctua- 
tions grow as described by the linearised fluid theory. 
Self-similarity applies to the full evolution to the extent 
that this self-similar temporal evolution at linear scales 
becomes "imprinted" on smaller non-linear scales. The 
mechanism by which this happens is simply the collapse 
of the initial mass fluctuations at large scales, at time 
scales fixed by linear theory. Self-similarity is thus a 
good approximation to the extent that the clustering am- 
plitudes at any scale depend only on the prior history of 
larger scales. In terms of power transfer in the evolution 
of clustering, self-similarity can thus be interpreted qual- 
itatively as indicating that there is a much more efficient 
transfer of power from large to small scales than in the 
opposite direction. Our results here show that this is true 
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FIG. 12: Collapse plot of A'^{k,t) for the SL128 simulation 
using as scaling factor Rs{t) as described in Eq.(^7[). 



also in more "blue" initial conditions with a small k PS 
P{k) cx fc" and n > 1, in which the variance of mass in 
real space spheres is dominated by fluctuations at much 
smaller scales which evolve in the non-linear regime. 

These points are further illustrated in Fig. |l^, which 
shows a "collapse plot" for the temporal evolution of 
A'^{k,t) = k^P{k). It follows from Eq. @) that, when 
self-similarity applies, we have the behavior 



A\k,t) = A^iR4t)k,ts) 



(37) 



where, as above, ts < t is an arbitrar y in itial time and 
Rsit) is given by Eq. (H). In Fig. || is plotted the 
rescaled function at each time, starting from tg = 0. At 
small k, we see that right from the initial time the self- 
similarity is indeed followed (as the rescaled curves are 
always superimposed at these scales) . This is simply be- 
cause linear theory, which is valid at these scales, gives 
such a behavior. As time progresses we see the range of k 
in which the curves are superimposed increases, extend- 
ing into the non-linear regime. Thus the self-similarity 
"propagates" progressively from small k to larger k, car- 
ried by the scales which are evolving non-linearly. Note 
that the behavior at asymptotically large k is simply 
A'^{k,t) cx k^/riQ (where no is the mean particle den- 
sity) at all times, corresponding to the shot noise present 
in all particle distributions with average density no and 
which by definition does not evolve in time^**. 

Let us finally return to the question of the breakdown 
of linear theory. In our discussion of Fig. in Sec. IIIB 1 



^"^ Note that in the two-point correlation function this time inde- 
pendent discrete contribution appears as a singularity at the ori- 
gin. It therefore does not "pollute" the collapse plots for ^(r, i). 
This also explains why one is able to identify the scaling behav- 
ior more readily by eye in this quantity. The collapse plot for 
a-^{r,t), which we have not shown, is similar to that for A^(fc). 
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FIG. 13: Behavior of A^{k,t) in SL128 together with the 
prediction of Unear theory (LT) . The points correspond to the 
value of A'^{kt,t), at times (from right to left) 1, 2, 3, 4, 5, 6 
and 7, where fc* is the wave-number above which the evolution 
the PS is no more longer well approximated by linear theory. 



above, we noted that the scale-independent amplification 
of linear theory describes very well the evolution of the 
PS up to a wavenumber k, which we denoted and 
which decreases with time. A question of interest is what 
the criterion is which determines this scale, i.e., what the 
criterion is for the application of linear theory. We can- 
not answer this question rigorously without considering, 
at least, the next order in this perturbative treatment 
We will not attempt to do so here, but rather consider de- 
termining such a criterion phcnomenologically (i.e. from 
the simulations). 

In principle this criterion may, in general, be quite 
complicated, as it would be expected to depend on the 
fluctuations present at all scales. Once we are in the 
self-similar regime, however, we expect that all charac- 
teristic scales in k space, and in particular k^(t), should 
scale oc Rj^{t). Such a time dependence results if one 
supposes kf{t) determined by a dimensionless quantity 
having some given amplitude. The evident simple crite- 
rion which suggests itself is 



A^(fc*(t),<) = constant. 



(38) 



Fig. shows the evolution of A^(fc,t), together with 
the evolution in linear theory. The points (small black 
circles) mark the approximate value of k^ at each time, 
determined as the scale at which the full evolution devi- 
ates from the linear theory in each case^^. The horizontal 



It is in fact possible [ [u| to write the equation for the evolution 
of density fluctuations in k space in a convenient form for this 
purpose, with all corrections to the linearised fluid limit in two 
formally simple terms. See also H] for discussion of these issues. 
The value of A;* in Fig. (used for the points) has been estimated 



line shows that, starting from about i = 3, when the self- 
similarity has set in, Eq. ( ^8|) with the constant set equal 
to unity is a reasonably good fit to the observed fc*(i). 
The deviation of the last point, at t = 7 can be attributed 
to finite size effects, as we see that at this time the small- 
est k modes in the box are no longer described well by 
the linear evolution. 

For n < 1, because of Eq. (p3|), the criterion Eq. ( |3^ ) 
for the breakdown of linear theory is equivalent to one 
stating it as a threshold value of the real space variance. 
In the current case, with n > 1, there is no such evident 
equivalence, as the mass variance at scales r > k~^(t) 
are not directly determined by the fluctuation of these 
modes, but rather by the fluctuations in larger k modes. 
Thus in this case the physical criterion for the breakdown 
of linear theory is really more appropriately given in k 
space 



6. Role of two-body correlations 

The gravitational force on a particle in an SL is dom- 
inated, for small S, by that exerted by its six nearest 
neighbors, and for large 6, by its single nearest neighbor 
[ psf . One thus expects that, at sufficiently early times, 
the dynamical evolution should be well approximated by 
neglecting all but these dominant contributions to the 
force. It has in fact been shown in that the early time 
evolution of simulations of small i5 SL can be well approx- 
imated by a two phase model: in a first phase the particle 
moves under the effect of its six nearest neighbors, and 
then subsequently, when the lattice symmetry is broken, 
under the effect only of a single nearest neighbor. The 
first non-linear correlations then emerge as these nearest 
neighbor pairs fall toward one another. 

As described in Sec.|l| the relation Eq. ( |l8| ) holds in the 
approximation that the correlations are primarily due to 
correlated pairs of nearest neighbor particles. Its valid- 
ity is thus a good probe of the probable adequacy of a 
dynamical model like that just described 
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In Fig 

shown the comparison of the two relevant quantities 
correlation function measured in the simulation and the 
one reconstructed by using Eq.(p^, i.e., by considering 
explicitly only nearest neighbors correlations. We see 
that the relation holds very well until t = 1. 

This is a quite striking and surprising result: the form 
of ^(r, t) — which is subsequently that which is observed 
to scale in the asymptotic self-similar evolution — has 



using the following criterion: | In P(A:, , t)/ In PltC^* , — 1| = 

0. 05, where PurikuS) is simply the initial PS amplified by linear 
theory, i.e., Eq. (|24[). 

If one wishes to define a real-space scale directly, this can be 
done by using the mass variance defined in a Gaussian window, 

1. e., with W{k;r) in Eq.nfn) given by a Gaussian of width 1/r. 
This is really just a trivialway of restating the k space condition 
in real space. 
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FIG. 14: Two-point correlation function in SL128 at times 
t — 0.5, 1, 1.5 (thin lines) together with the approximation 
got from the NN PDF (thick lines). For clarity the behaviors 
at different times have been arbitrarily rescaled on the i-axis. 



already emerged at a time when nearest neighbor inter- 
actions play a crucial role in the dynamics. In the previ- 
ous sections, however, we have seen that this asymptotic 
behavior is characterised by a time dependence derived 
in a fluid hmit. Such a limit is normally expected to be 
valid in the opposite case that two or few body inter- 
action with nearest neighbor particles can be neglected 
(rather than being dominant in the approximation just 
considered). Indeed we noted that when the asymptotic 
scaling behavior there is necessarily no explicit depen- 
dence on the characteristic length scales in the system 
— and notably those associated with the discreteness of 
the distribution which directly enter in determining the 
strength of nearest neignbour forces. We will discuss this 
point further below after a presentation of results of the 
other SL simulations we have performed. 



7. Dependence on the normalized shuffling parameter 

As discussed above SL initial conditions may be char- 
acterized, for their gravitational evolution, by the single 
dimensionless parameter S. Our analysis until now has 
concerned solely the simulation SL128 and thus only a 
single value {S = 1). Our primary result — that this 
system tends in a few dynamical times to a "self-similar" 
evolution — would be expected to be true for any (fi- 
nite) value of 6: this particular spatio-temporal scaling 
behavior is determined solely by the fc^ form of the small 
k PS, which is invariant under changes in i5. Thus the 
only thing that we would expect to change non-trivially 
when S changes is the transient regime to the asymptotic 
self-similar behavior. Specifically we might expect both 
the duration of this transient and its nature to change. 

The emergence of self-similarity in the evolution cor- 



responds, as we have discussed at length, to a behavior 
which is explicitly independent of the discreteness scale 
i characterizing its particle-Hke nature. The simplest in- 
terpretation of this behavior — and the usual one in cos- 
mology — is that this corresponds to a fluid-like behavior 
of the system i.e. to an evolution which can be described, 
at both linear and non- linear scales, by a set of non-linear 
fluid equations approximating the particle dynamics 
If this interpretation is correct, any (5-dependent effects 
in the evolution of SL with different (5, but with the same 
large scale fluctuations (i.e. small k PS) can then be con- 
sidered as "discreteness effects" . 

This equivalence of the fluid limit of the evolution from 
SL with different S can be seen even more explicitly as 
follows, for the case that S is small. In this limit the 
so-called Zeldovich approximation to the fluid limit evo- 
lution (see Appendix y| below) is valid. Each element of 
the fiuid then moves according to 

x(q,i) =q + /(i)u(q,t = 0) (39) 

where q is a Lagrangian (time-independent) coordinate, 
which we can take here to be the lattice point from which 
the particle/fluid element is displaced, and u(q, t — 0) is 
the displacement of the particle/fluid element at the ini- 
tial time. The function f(t) is simply the growth factor 
of the fluctuations in linear theory. The effect of the evo- 
lution, in this limit, is thus manifestly to transform one 
SL into another one with a different (larger) S. Thus in 
the linear fluid limit, starting from a small S, the evolu- 
tion of the system should be identical (statistically, and 
up to an overall scale transformation) to that of an SL 
with a larger 5. 

The simulations SL64, SL32, SL24, and SL16, as we 
have deflned them allow us to explore the S dependence 
(and thus non-fluid effects) in the evolution from SL ini- 
tial conditions. As described above in Sec.||, we have 
chosen in each case a combination of S and £ which leaves 
the amplitude of the PS constant (in the length units we 
have chosen, fixed by the box size in these simulations). 
This choice means that the evolution of any two simu- 
lations in time should agree (without any rescaling) if 
the evolution of both may be well described by the fiuid 
limit: this is governed, as we have seen, by the evolution 
of the fiuctuations at large scales which are identical. 

These statements are of course true in the approxima- 
tion that effects introduced by the finite box-size of the 
simulation, the softening of the force and any other ef- 
fects of the numerical discretization of the evolution, are 
negligible. We noted that in a finite box, and with soft- 
ened gravity, one has two additional parameters, which 
one can choose as £/L and e/L. The simulations SL64, 
SL32, SL24 and SL16 in Table | correspond, as we have 
described, to chosen fixed values of these two parameters. 



More precisely the system is assumed then to evolve as described 
by a set of Vlasov-Poisson equations. See appendices below. 
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In order to control for dependence on the box size L, we 
have chosen SL64 to have the same 6 as SL128, so that 
the two sets of initial conditions differ only in the box 
size. Thus these two simulations should give precisely 
the same (averaged) results as long as finite size effects 
play no significant role. We will not report in this paper 
the sensitivity of results to the choice of e. We have, how- 
ever, verified that, for a considerable range of variation of 
£ to smaller values than the one used in the simulations 
we report here, there is no notable effect on our results 
in the range of scales r > e where we assume them to be 
valid. 

In Fig. |l5| are shown ^(r, t) in each of the five SL simu- 
lations given in Table |, for different times, starting from 
t = 1 when the structures first develop in the largest sim- 
ulations until t = 6 when the scale of homogeneity has 
reached a significant fraction of the box size. 

The first point to note is the excellent agreement be- 
tween the results of SL128 and SL64, which differ only in 
the size of the simulation box. This assures us that the 
finite size effects due to the different values of i, up to 
the time we have shown, are very small in the SL64 simu- 
lation, and we will assume the same is true for the SL32, 
SL24, SL16 simulations (in ascribing the differences be- 
tween them solely to the change in S and not to that in 
the number of particles). 

Considering now the evolution of the other simulations 
we observe that: 

• As ^ decreases the time increases at which the sys- 
tem begins to evolve and form strong non-linear 
correlations (i.e. develop a region with ^(r, t) ^ 1). 
This is a qualitative behavior expected also in the 
fluid limit: in the Zcldovich approximation Eq.(|39|) 
the displacements grow at a rate given by the func- 
tion f{t) which is independant of scale. Thus, 
starting from a smaller S, the time at which non- 
linear structures form (when (5 ^ 1) is necessarily 
longer^^. 

• When the first non-linear correlations develop there 
is a manifest d dependence in the correlation func- 
tions, i.e., the correlations are not (statistically) 
equivalent to those in the larger S simulations. This 
means that at the time these correlations emerge, 
the smaller 6 system is not evolving as in its fluid 
limit. If it were it would be in agreement with the 
larger 5 simulation with the same initial power at 
the relevant scales. 

• Initially the non-linear correlations formed in each 
system "lag behind" those in the larger 6 simula- 
tion, i.e., ^(r, t) typically has a smaller amplitude in 



Equivalently one can say that the larger I system is "missing in- 
put power" above its Nyquist frequency compared to the smaller 
I simulation. 



the smaller 6 simulations. As it evolves the smaller 
S system eventually "catches up" with the larger 
ones, its correlations eventually agreeing very well 
with those in all the larger 6 systems over a signif- 
icant range of scale. 

• The form of the non-linear correlation function in 
the asymptotic regime — the self-similar regime we 
have discussed above — emerges to a very good 
approximation at a time when there is still a quite 
visible "lag" in amplitude. 

In Fig. |l^we show also the evolution of i?s(i), inferred 
in each case, as in our analysis of SL128 above, by the 
determination of the factor which describes the spatio- 
temporal scaling once it emerges as a good approxima- 
tion. We observe that in each case we have, as for SL128, 
a regime of approximate spatio-temporal scaling of the 
non-linear correlation function before the asymptotic self- 
similar regime is reached. In this regime Rs{t) is smaller 
in amplitude than in the asymptotic regime, correspond- 
ing to the "lag" of the smaller 6 simulations described. 
However Rs{t) evolves more rapidly than in the asymp- 
totic regime, allowing each system to "catch up" with the 
(^-independent behavior. Note that these observations 
again confirm that the corresponding regime in Rs{t) in 
SL128 should indeed not be ascribed to a first self-similar 
phase driven by the Poisson fiuctuations present in this 
case. 

Both this "lagging" and the role of nearest neighbor 
interactions in the formation of the first structures can 
be explained in the framework of a refined version of the 
"two phase model" of j2^ for the early time correlations. 
We will present this model in detail elsewhere, and re- 
strict ourselves here to a few qualitative comments. 

A very good approximation to the evolution of a per- 
turbed lattice is provided by a perturbative treatment 
described in Q. The force acting on particles is 
written as an expansion in the relative displacement of 
particles, in a manner completely analogous to a stan- 
dard technique used in solid state physics to treat per- 
turbations to crystals. One can then do a linear mode 
analysis in k space to determine the eigenmodes of the 
displacement fields under gravity. While at small k one 
recovers the simple k independent amplification of linear 
fluid theory, the effect at larger k (i.e. k ~ fc^v) is, for all 
but some very specific modes, to slow down the growth 
of fiuctuations. Thus the "collapse time" for fluctuations 
at scales of order the inter-particle distance are indeed 
slowed down, as observed here, compared to linear fluid 
theory. 

This approximation to the early time evolution breaks 
down when the force on a particle starts to be dominated 
by a single nearest neighbor. At this point particles start 
to accelerate toward their neighbor, giving rise to strong 
two-body correlations which are, as we have seen above, 
the dominant contribution to the measured two-point 
correlations at non-linear scales at early times. We have 
remarked that, given this manifestly "non- fluid" mecha- 



18 




10" 

lo' 

10° 
10-' 
10-^ 
10-3 
10-^ 







SL128 






SL64 






SL32 






SL24 






SL16 








e 


1 111 , \ 







FIG. 15: Evolution of the two-point correlation function in the different simulations at times 1, 2, 3, 4, 5 and 6. The four thick 
arrows represent the different mesh sizes t while the thin one corresponds to the value of the softening length e. 



nism for the formation of these correlations, it is some- 
what surprising to see approximately the same two-point 
correlations maintained in the "self-similar" regime, if 
this regime is interpreted as the result of a purely fluid- 
like evolution. Two possible, but very different, explana- 
tions for this are the following: 

• The self-similar evolution of the system in the non- 
linear regime is not correctly interpreted as a man- 
ifestation of a purely fluid limit of the N body sys- 
tem. Its timcscalcs arc dictated by the fluid limit 
(giving the collapse time for fluctuations at large 



scales), but its non- linear dynamics are intrinsically 
discrete; 

• The early time non-linear correlations, well de- 
scribed by a discrete dynamics, approximate well 
those in the fluid limit because the non-linear fluid 
dynamics is in fact physically well approximated 
by a discrete system, i.e., the non- linear evolution 
of the fluid, in the relevant phase of moderately 
strong non-linear correlations (^(r) < 10^), is well 
described as the evolution of "lumps" of fluid to- 
ward "nearest neighbor lumps". 
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FIG. 16: Evolution of the rescaling factor Rs{t) in the dif- 
ferent simulations. Also shown is the self-similar behavior 
Eq. (H). 

We will evaluate these two quite different interpreta- 
tions of our results more quantitatively in future work. 



IV. DISCUSSION AND CONCLUSIONS 

To conclude we first summarize our conclusions, and 
then make a few remarks on open questions to be ex- 
plored in further works. 

We have studied the evolution under their Newtonian 
self-gravity, in a static euclidean space, of classical point 
particles initially distributed in infinite space in a quasi- 
uniform manner. This is a paradigmatic problem of the 
out of equilibrium statistical mechanics of long range in- 
teracting systems, which has received little attention in 
this context. Specifically we have considered a one rel- 
evant parameter class of initial conditions in which the 
particles are randomly perturbed off a lattice. We have 
found that our simulations converge aysmptotically (but 
for times smaller than those at which the size of the fi- 
nite simulation box becomes relevant) to solutions char- 
acterized by a simple spatio-temporal scaling relation in 
which the temporal dependence of the scaling can be de- 
rived from the linearized fluid theory. These results are 
qualitatively very similar to those observed in numerical 
studies in the context of cosmology, i.e., for expanding 
space-times and for more complex initial conditions in 
which the displacements of the particles off the lattice 
are correlated in order to produce the PS of fluctuations 
of cosmological models. More specifically, the observed 
spatio-temporal scaling is a simple generalization of what 
is known in the cosmology literature as "self-similarity" 
in an expanding universe to the case of (i) a static uni- 
verse, and (ii) a PS P{k) oc k'^. Further we have ob- 
served that there is a transient phase to this behavior, in 
which already, to a good approximation, the same spatio- 
temporal scaling relation holds for the two-point correla- 



tion function ^(r, t), but with a more rapid temporal evo- 
lution of the scaling factor. We have noted that the lag- 
ging of the evolution behind the asymptotic behavior in 
this regime can be ascribed to effects of discreteness (i.e. 
non fluid effects) slowing down the evolution of fluctua- 
tions at scales comparable to the inter-particle distance 
which have been quantified in ^ . We have seen also 
that the form of the correlation function emerges already 
at the very early times when the first non-linear correla- 
tions develop due to two-body correlations which develop 
under the effect of nearest neighbor interactions. 

The gravitational evolution of a SL in a static universe 
thus shares the qualitative features of similar, but more 
complicated models, in cosmology. It thus provides a sim- 
plified "toy model" in which to study some fundamental 
problems which remain open concerning the evolution of 
these systems, which have been studied extensively in 
numerical simulations but remain poorly understood an- 
alytically, notably: 

• The absence of a theory which adequately explains 
the shape (i.e. functional form) and evolution of 
the observed non-linear correlations. 

• The absence of a "theory of discreteness errors". 
In cosmology simulations of particles displaced off 
lattices (or "glasses") aim to reproduce the evolu- 
tion of a self-gravitating fluid. There is currently 
very little systematic understanding of how well 
this evolution is actually approximated. We have 
highlighted in this paper that the SL gives a very 
well defined, and simplified, framework in which to 
address this problem. 

Let us remark finally on a few other points: 

• We have worked here with initial velocities set equal 
to zero. In exploring the analogy with cosmological 
simulations there is another choice of initial veloc- 
ities which is natural. This is that corresponding 
to that given by the Zeldovich approximation dis- 
cussed above, with f{t) chosen in Eq. ( p^ to corre- 
spond to the purely growing mode of density fluc- 
tuations, i.e., f{t) — e*/'^''y". The initial velocities 
are then simply the initial displacements divided 
by Tdyn. This introduces no further new charac- 
teristic scales in the initial conditions. Its effect 
on the evolution will be to make the transient to 
self-similarity slightly shorter, but it will not signif- 
icantly change any of our flndings or conclusions. 

• We have made a specific choice of PDF for our shuf- 
fling, given in Eq. (^. We expect different choices 
again to modify slightly the nature of the transient, 
but not the self-similarity. This latter, as we have 
emphasized, depends only on the form of the 
PS at small k, which is in fact the same for any 
PDF with finite variance. Indeed the coefficient 
of the k'^ is just given by this variance, and the 
difference between PDFs will manifest themselves 
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in modifications of the fluctuations at small scales 
(i.e. larger k). For example if the two PDF have 
different fourth moments, this will be reflected in a 
different coefficient in the k'^ correction to the small 
k PS. Just as in the case of velocities, there is a nat- 
ural choice if one wishes to maximize the analogy 
with cosmological simulations: a simple Gaussian 
PDF which is what is used in this context. In fact 
this choice is also natural from another point of 
view, as we will explain in detail in a forthcoming 
article p5|: when one considers constructing new 
particle distributions by a simple "coarse-graining" 
on some scale, the SL with Gaussian PDF, due to 
the Gentral Limit Theorem, has the property of be- 
ing the unique one which is invariant under such a 
coarse-graining. 

• We have reported in this paper simulations in which 
the softening e has been kept fixed (in our chosen 
length units). We have mentioned that we have 
checked that our results for clustering amplitudes 
above this scale are robust to the use of significantly 
smaller values of e. A more extensive and system- 
atic study of the role of this parameter would, how- 
ever, be of interest, specifically with the goal of un- 
derstanding in detail how the clustering properties 
are modified by it at small scales. 
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sure. The set of equations closes if p(x, t) is specified as 
a function of the density. 

As it is shown in |l^ , this set of equations can be 
obtained after certain approximations from the Vlasov 
equation coupled to the Poisson equation: 

[a* + V • Vx - Vv$ • Vv] /(x, V, t) = (A2a) 

where $ satisfies the modified Poisson equation 



V2$(x,<) = AttG 



I I /(x,v,t)d3v 

JR3 



PO 



(A2b) 



and f{ii.,\f,t) is the density of particles in the infinitesi- 
mal volume d'^xd^v at (x, v) at time t. These equations 
can themselves be derived as truncations of a BBGKY 
hierarchy p^, pj[ |, or starting from a Liouville equation 
for the full ("spiky") one particle phase space density 

gig. 

By performing a perturbation analysis ,for the case of 
pressureless (i.e. highly non-relativistic or "cold") mat- 
ter, around p = po and v = with the set of Eq. (|Al|), 
one finds at first order that the evolution of the density 
contrast (5(x, t) = (p — po)/po is described by the differ- 
ential equation 



S{x,t) = 47rGpo(5(x,t) 



(A3) 



or equivalently that each Fourier mode (5(k, t) evolves in- 
dependently of the others: 

^(k, t) = 4TTGpoS{k, t) . (A4) 
The general solution of Eq. (^) is §13] 

(5(x,t) = A(x) (^s/AttGpo + B(x) exp [-^AnGpo t) • 

(A5) 

For the case, as in this paper, in which the initial velocity 
is set equal to zero, one obtains 

(5(x, t) = <5(x, 0) cosh {^^AttGpq . (A6) 



APPENDIX A: FLUID EQUATIONS AND FLUID 
LINEAR THEORY 



The equations which describe the evolution of a sclf- 
^ra vitating fluid are the following (e.g. chap. II] or 
E4, chap. 5.2]) 



dtp + Vx • (pv) = 
(9tv + (v • Vx)v = g VxP 



Vx g = 
Vx X g 



~4ttG{p~ Po) , 
, 



(Ala) 

(Alb) 

(Ale) 
(Aid) 



where p(x, t) is the mass density, v(x, t) the velocity 
field, g(x, i) the gravitational field and p(x, t) the pres- 



APPENDIX B: LAGRANGIAN FLUID THEORY 
AND THE ZELDOVICH APPROXIMATION 

The previous appendix uses the Eulerian formalism of 
fiuid mechanics, in which one describes the evolution of 
the different quantities characterizing the fluid (veloc- 
ity, density and pressure) at each point of a flxed ref- 
erence frame. In the alternative Lagrangian formalism 
(see, e.g.,r 
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50|), one describes the evolution of the 
fluid in terms of the displacements of its elements with 
respec t to a reference frame. In that case, the equa- 
tions (Al) which describes the density and the velocity 
are then transformed into a set of equation describing the 
evolution of a displacement fleld f(X,t). The position x 
of a fluid element at time t is then written as 



x = f(X,t) 



(Bl) 
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where the coordinate X labels the fluid element consid- 
ered. One can choose this coordinate as the position of 
the fluid element at the initial time (which we assume to 
be 0): X = f (X, 0). The equations for f (X, t) in the case 
of a gravitating fluid can be found in, e.g., Note 
that in this reference, a fluid in an expanding universe is 
considered. The static case can be recovered by setting 
the expansion factor a{t) = 1 at all times. 

As in the Eulerian approach, one can perform a per- 
turbation theory in the Lagrangian approach. One writes 
f (X, t) ^ X + p(X, t) with p(X, 0) = 0, and performs 
a Taylor expansion in powers of p. At linear order in 
p(X,t), one obtains the following set of equations: 



V-(p- 
V X p 



AnGpop) 
-- 



-47rGpo(5(X, 0) 



(B2) 
(B3) 



where (5(X, 0) is the density contrast at i = 0. Writing 
the vector field p as the sum of a curl-free part pu and 
a divergence-less part pn (i.e. p^ can be written as the 
gradient of a scalar function, and pn as the curl of a 
vector field), one finds after some calculation that 



p(X,t) =p(X,0) 



Pd(X,0) 



[cosh (V47rG>o t) - 1] 
sinh {^/AirGpo t) 



Pii(X,0)t 

m- 



with the initial condition p(X,0) = 0. Note that 
p(X,0) = pd(X,0) since the gravitational force is con- 
servative. These Eqs. (B4) correspond to Eqs. (6) and 
(7) in m, with a{t) 1 and A = AnGp^. 



The asymptotic behavior of the solution (B4) is 



p(X,<) 



p(X,0) , Pi3(X,0) 



47rG'/9o V47rGpo 



exp \ \J AnGpo t 
(B5) 

By choosing pfl(X,0) = and p(X, 0)\/47rGpo = 
p(X,0), the solution is then directly in its asymptotic 
regime. This is the static space equivalent of the Zel- 
dovich approximation in an expanding background 

The linear approximation of the Lagrangian approach, 
which leads to the Zeldovich approximation as we have 
described, has proven to be very useful in the problem 
of gravitational clustering. With respect to the linear 
Eulerian approach, it has the advantage that it can de- 
scribe the evolution of density fiuctuations with a density 
contrast much greater than unity. 
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